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1.  Introduction 


Recent  interest  in  high  Mach  number  flow  has  pointed  out  deficiencies  of  classical  central  difference 
schemes  near  strong  discontinuities.  While  these  schemes  perform  admirably  in  sub-,  trans- 
and  even  low  supersonic  flows  [1,  2,  3,  4,  5],  problems  arise  due  to  oscillations  which  may  grow 
unbounded  in  the  vicinity  of  strong  shocks.  Inaccuracies  due  to  the  excess  artificial  dissipation 
required  to  suppress  this  behavior  led  to  the  development  of  a  variety  of  upwind-biased  schemes. 
These  schemes  rely  on  local  wave  propagation  theory  in  differencing  the  convective  terms  of  the 
governing  equations  throughout  the  domain. 

The  focus  of  the  current  effort  is  to  investigate  the  characteristics  of  several  upwind  schemes 
for  problems  of  realistic  complexity  in  high  speed  flows  against  the  backdrop  of  classical  central 
difference  techniques.  Upwind  schemes  for  the  Euler  and  Navier-Stokes  equations  may  be  distin¬ 
guished  by  the  approximate  Riemann  solver  used  to  evaluate  the  fluxes  at  the  cell  faces  and  the 
mechanism  through  which  they  obtain  higher  order  spatial  accuracy  [6,  7,  8]. 

The  approximate  Riemann  schemes  of  Roe  [9],  Osher  [10],  van  Leer  [11]  and  Steger- Warming  [12] 
currently  enjoy  wide  use  within  the  research  community.  These  solvers  are  extended  to  higher 
than  first  order  through  either  a  MUSCL  [13]  or  non-MUSCL  approach  [8,  14].  For  a  non-linear 
scalar  equation  or  a  linear  system  of  equations,  schemes  developed  with  either  approach  may  also 
be  TVD  through  the  incorporation  of  an  appropriately  constructed  flux  limiter. 

The  concurrent  development  of  many  of  these  algorithms  has  led  to  several  comparisons  of 
similar  schemes  [13,  15,  16]  or  comparisons  of  a  variety  of  schemes  for  simple  problems  [17, 18,  19, 
20,  21,  22].  However,  little  effort  has  been  made  to  study  a  wide  class  of  schemes  in  a  controlled 
setting  on  rigorous  problems. 

This  led  to  the  current  work  which  examines  several  MUSCL  and  non-MUSCL  schemes  with 
a  variety  of  Riemann  solvers  (Table  1.1)  for  a  realistically  complex  configuration.  These  schemes 
are  compared  with  central  schemes  using  second  and  fourth  order  artificial  damping.  The  shock 
on  cowl  lip  problem,  which  was  investigated  experimentally  by  Holden  et  al.  [23],  has  received 
much  attention  in  recent  literature  [24,  25,  26,  27]  and  provides  a  rigorous  test  case  of  practical 


1 


importance.  Edney  [28J  classified  shock  interference  patterns  into  six  types.  The  actual  pattern 
established  in  any  particular  case  depends  primarily  upon  the  strength  and  location  of  the  im¬ 
pinging  shock,  the  characteristics  of  the  incoming  flow  and  the  cowl  shape.  Edney  predicted  the 
highest  mechanical  and  thermal  loading  in  the  Type  IV  interaction  which  establishes  itself  when 
the  impinging  shock  intersects  the  subsonic  portion  of  the  cowl  bow  shock.  Figure  1  contains 
a  schematic  of  the  Type  IV  interaction  that  was  investigated  experimentally  by  Wieting  and 
Holden  [29].  A  supersonic  viscous  shear  layer  emanates  from  the  point  of  impingement.  This 
layer  forms  a  “jet”  which  strikes  the  body  after  passing  through  a  terminating  bow  shock.  Since 
inviscid  features  dominate  this  flow,  the  investigation  is  free  to  focus  on  the  discretization  of  the 
inviscid  terms  only.  This  problem  includes  strong  shocks  and  shear  layers  interacting  in  the  pres¬ 
ence  of  a  body.  These  phenomena  create  many  length  scales  in  the  flow  and  the  shocks  and  shear 
layers  are  not  necessarily  aligned  with  the  grid  [23].  Furthermore,  under  certain  circumstances,  a 
weak  instability  arises  in  the  interaction  which  may  or  may  not  be  damped  out  by  the  truncation 
error  of  a  given  numerical  method  [25,  30]. 

The  present  work  compares  all  schemes  on  identical  grid  systems  at  three  levels  of  refinement. 
Particular  emphasis  is  placed  on  convergence  and  accuracy  of  the  algorithms  under  consideration. 
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Scheme  No. 

Spatial  Discretization 

1 

Central  differencing  (5-point  operator)  +  2nd  and  4th  order  dissipation 

2 

Central  differencing  (9-point  operator)  +  2nd  and  4th  order  dissipation 

3 

Steger- Warming  flux  vector  splitting  +  MUSCL 

4 

van  Leer  flux  vector  splitting  +  MUSCL 

5 

Roe  flux  difference  splitting  +MUSCL 

6 

“Symmetric”  TVD  of  Yee  (non-MUSCL) 

7 

“Upwind”  TVD  based  on  Harten  and  Yee  (non-MUSCL) 

Scheme 

Limiter 

Time  Integration  Scheme 

Reference 

1 

— 

Explicit  Runge-Kutta  scheme 

[31] 

2 

— 

Explicit  Lax  Wendroff 

[32] 

3 

Minmod  limiter 

Gauss-Seidel  line  relaxation 

[12,  33] 

4 

Minmod  limiter 

Gauss-Seidel  line  relaxation 

[11,  33] 

5 

van  Albada  limiter 

Explicit  Runge-Kutta  scheme 

[9,  8] 

6 

Limiter  (Ref.  [8],Eqn.  3.51h) 

Explicit  Runge-Kutta  scheme 

[34,  35] 

7 

Limiter  (Ref.  [8],Eqn.  3.51f) 

Explicit  Runge-Kutta  scheme 

[36,  35] 

Table  1.1:  Brief  description  of  schemes  and  their  nomenclature 
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2.  Brief  Description  of  Schemes  Investigated 


In  descretizing  the  Navier-Stokes  equations,  virtually  all  schemes  rely  on  a  centered  approximation 
of  the  viscous  terms.  In  an  effort  to  focus  solely  on  the  convective  flux  approximation  of  these 
upwind  schemes,  this  investigation  considers  the  Euler  equations.  Additionally,  excluding  the 
viscous  terms  removes  the  diffusive  loss  mechanism  from  the  governing  equations  and  provides  a 
homenthalpic  reference  state  which  may  be  used  to  assess  the  accuracy  of  the  various  methods. 
Writing  the  Euler  equations  in  conservation  form: 


dF  dG 
dx  dy 


=  0 


(2.1) 


U  is  the  solution  vector  of  conserved  variables  and  F  and  G  are  the  flux  vectors.  The  gas  is 
assumed  to  be  calorically  perfect  governed  by  the  equation  of  state  p  =  pRT  where  p  is  the 
pressure,  p  the  density,  R  the  gas  constant  and  T  the  static  temperature. 

This  study  focuses  on  the  performance  of  MUSCL  and  non-MUSCL  higher  order  upwind 
schemes.  Suppressing  the  vector  notation,  denoting  UL  and  UR  as  the  left  and  right  states,  a 
brief  description  of  each  algorithm  follows  with  reference  to  the  evaluation  of  the  inviscid  flux  at 
the  cell  surface  i  + 


MUSCL  Type  Upwind  Schemes 


The  MUSCL  type  schemes  considered  include  Steger- Warming  and  van  Leer  flux  vector  splitting 
and  Roe  flux  difference  splitting. 


•  Steger- Warming  flux  vector  splitting  [12]: 

The  Steger- Warming  flux  vector  splitting  scheme  may  be  written  as: 

Fi+i  =  F+(t£j)  +  F-{V*  j)  (2.2) 

With  A  =  OF/dU, 

A  =  R~lAR  =  R~l(A+ +  A~)R  =  A+  +  A~  (2.3) 
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F+(U)  =  A+U,  F~(U)  =  A~U  (2.4) 

Here  A  is  a  diagonal  matrix  consisting  of  the  eigenvalues  of  A,  A+  and  A-  denote  the  splitting 
of  the  eigenvalues  into  positive  and  negative  components  and  R  is  the  similarity  transformation 
diagonalizing  A,  The  MUSCL  approach  [13]  brings  UL  and  UR  to  second  order  accuracy: 

(2'5> 

«£}  =  +  (2.6) 

where: 

Ai+i  =  minmod(A<+i,At_i)  (2.7) 

and  A  - ,  i  =  UR  -  UL.  With  the  addition  of  the  minmod  limiter,  the  algorithm  reverts  to  first 
*+3 

order  accuracy  at  shocks  in  order  to  preserve  monotonicity  within  of  the  solution. 


•  van  Leer  flux  vector  splitting  [11]: 


The  functional  form  of  the  van  Leer  scheme  is  identical  to  the  Steger- Warming  algorithm 
(Eqn.  2.2)  and,  in  fact,  for  supersonic  flow  the  evaluation  remains  unchanged.  For  subsonic  flow, 
were  revised  by  van  Leer  to  avoid  the  discontinuity  exhibited  across  sonic  lines  by  the  Steger- 
Warming  algorithm.  Denoting  the  Mach  number  normal  to  the  interface  as  Mx,  for  subsonic  flow 
(|M,|  <  1), 


F±  = 


±pa[\{Mx  ±  1)}  =  ft 

f-r-l)ti±2a 

/,*» 

if  + i\ 


(2.8) 


This  scheme  also  obtains  second  order  accuracy  with  the  MUSCL  approach  (Eqns.  2.5  through  2.7). 
Both  Steger- Warming  and  van  Leer  methods  utilize  the  cell-centered  finite  volume  approach  in 
transformed  space  [37]. 
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•  Roe's  flux  difference  splitting  [9]: 


The  Roe  flux  difference  splitting  reads  as  follows: 


(2.9) 


where  ( * )  indicates  evaluation  at  the  Roe  averaged  state  between  UL  and  This  scheme 

extends  to  second  order  accuracy  through  the  MUSCL  approach  (Eqns.  2.5  and  2.6)  and  was 
applied  in  conjunction  with  the  van  Albada  limiter  [38,  8]. 


A. 


*♦»- 


-} 


(A,+  {)"  +*!l  +  Ai+}  [ (A.'-})3  +  *» 


ALi  +  A?-i.  +  2«i 

*  '  2  *  2 


(2.10) 


where  is  a  small  constant  to  prevent  division  by  zero.  In  contrast  to  the  above  flux  vector 
split  algorithms.  Roe’s  scheme  may  violate  the  entropy  condition  when  the  eigenvalues  at  the  Roe 
averaged  state  vanish.  Following  Harten  [39]  the  eigenvalues  |A|  of  |A|  are  modified  when  they 
are  below  some  small  threshold  8: 


"-VP 


|A|<« 


(2.11) 


Non-MUSCL  Approach 

The  works  of  Yee  [34]  and  Harten  and  Yee  [36]  provide  two  popular  examples  of  non-MUSCL 
upwind  schemes.  The  first  order  formulations  of  both  reduce  identically  to  the  Roe  solver. 


•  Yee’s  “Symmetric”  TVD  scheme  [34,  35]: 


This  method  takes  the  following  form: 

F.+l  =  +  f(K)l  +  j8j+i*i+i 


(2.12) 
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where  is  the  Ith  component  of  $  given  by: 

«+3 

[ai+$  ~  $!+ ■  J 

Here  again  |A|  is  thresholded  as  in  the  Roe  scheme  (Eqn.  2.11).  Also: 


(2.13) 


«!+J  =  A'+'i  W+.  -  Vi)  (2.14) 

is  essentially  the  difference  of  the  characteristic  variables  and  Ql  ,i  is  the  limiter  taken  from  [8] 
(Eqn  3.51h). 

=  minmod  [2a'._i,2a|+i,  ^  +  a'+i)]  (2.15) 

It  is  the  addition  of  this  limiter  which  brings  the  scheme  to  second  order.  Setting  Ql  =  0  reduces 
this  method  identically  to  Roe’s  first  order  scheme. 

•  Harten  and  Yee’s  “Upwind”  TVD  Scheme  [36,  35]: 

The  modified  flux  approach  of  Harten  and  the  modifications  by  Harten  and  Yee  provide  an 
alternative  second  order  approach.  This  may  be  couched  in  the  form  of  Eqn.  2.12  by  redefining 
$  for  the  fth  component. 

(*!+*)  (ff'+i + *9  -  *  (*!+ * + 4 *)  «!+*  (2-16) 


where: 


* 

- 1*  (*4)  ■ 


(2.16) 

if  a'+i  #  0 

if  a‘  =  0 

(2.17) 

and  enforces  the  entropy  condition  by  preventing  |zj  from  reaching  zero. 

fW  W>* 


M 

25 


\z\  <6 


(2.18) 


The  limiter  function  g  is  taken  from  [8]  (Eqn.  3.51f). 


9'  a!_i + ai+i + ^ 


(2.19) 
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The  three  previous  methods  (Roe  flux  difference  splitting.  Symmetric  and  Upwind  TVD)  are 
implemented  in  a  cell  vertex  flnite  volume  formulation  [40]. 

Central  Difference  Schemes 

To  provide  a  datum  by  which  to  evaluate  the  upwind  procedures,  the  comparison  includes 
formulations  of  two  popular  centered  schemes. 

•  Jameson’s  classical  five-point  scheme  [31]  with  second  and  fourth  order  artificial  dissipation 
implemented  in  a  node  based  manner  [40]. 

•  Ni’s  nine-point  finite  volume  formulation  of  the  classical  Lax-Wendroff  procedure  [32]  with 
the  addition  of  second  and  fourth  order  dissipation  [41]. 

Three  distinct  time  integration  techniques  are  utilized  in  this  research.  The  Steger- Warming  and 
van  Leer  flux  split  schemes  are  integrated  in  time  using  the  implicit  Gauss-Seidel  line  relaxation 
algorithm  described  by  MacCormack  [33].  Following  the  conclusions  of  Liou  and  van  Leer  [18],  the 
implicit  operator  utilizes  the  Steger- Warming  Jacobians  for  both  splittings.  With  the  exception 
of  the  Lax-Wendroff  scheme,  all  other  schemes  use  explicit  multistage  Runge-Kutta  algorithms. 
Table  1.1  highlights  the  schemes  and  numbers  them  for  discussion  in  later  sections. 
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3.  Overview  and  Classification  of  Computations 


In  order  to  assess  the  quality  of  the  current  implementation  of  each  scheme  for  high  speed  flow, 
some  basic  issues  are  examined  on  a  simple  problem.  A  noninterfering  2-D  cylinder  (Case  A) 
at  Mach  8.03  provided  a  qualitative  check  of  the  convergence  properties  and  the  shock  stand  off 
distances  obtained  by  each  scheme.  This  case  also  brought  to  light  the  behavior  of  the  sonic  line 
in  view  of  the  limiting  characteristic  theory  provided  in  classical  references  [42]. 

This  simple  problem  gives  some  preliminary  insight  into  the  shock  capturing  capability  and 
overall  accuracy  of  the  various  schemes.  The  shock  capturing  capability  is  investigated  by  line  plots 
across  grid-aligned  and  nonaligned  portions  of  the  bow  shock.  Shock  resolution  is  also  examined 
as  a  function  of  Mach  number  with  the  addition  of  a  Mach  16.33  noninterfering  cylinder.  In 
addition,  a  grid  refinement  study  on  stagnation  enthalpy  losses  along  a  streamline  at  Mach  8.03 
provides  an  indication  of  the  net  error  in  each  scheme. 

With  results  of  these  basic  investigations  in  mind,  the  discussion  turns  to  the  complex  shock- 
shock  interaction  mentioned  earlier.  This  problem  (Case  B)  results  in  severe  mechanical  and 
thermal  gradients  throughout  the  field  and  on  the  boundary.  In  the  sketch  of  Figure  1,  the 
terminating  normal  shock  results  in  surface  pressures  which  exceed  the  non-impinging  case  by  an 
order  of  magnitude.  Numerically,  this  case  focuses  attention  on  scheme  performance  in  complex 
flows.  Moreover,  since  the  geometry  remains  simple,  this  complexity  stems  from  shock  and  slip 
line  structure  within  the  flow  field,  rather  than  from  an  arbitrary  boundary  geometry.  This 
structure  gives  rise  to  many  interdependent  time  scales  which  require  many  characteristic  times 
to  pass  before  a  steady  state  establishes  itself.  In  an  effort  to  add  rigor  to  this  comparison,  the 
Type  IV  problem  is  investigated  on  a  series  of  three  successively  refined  meshes.  All  flow  and  grid 
parameters  are  listed  in  Table  3.1. 
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Case 

Mach  Number 

Grids  ( IL  x  JL) 

Flow 

n 

8.03,16.33 

20  X  20,40  x  40 

Noninterfering 

B 

8.03 

60  x  40,120  X  80,240  x  160 

Interfering  (Type  IV) 

Table  3.1:  Classification  of  computations  for  flow  past  a  2-D  cylinder 
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3.1  Criteria  for  Comparison 


Throughout  this  effort,  emphasis  was  placed  on  proper  comparison  of  the  various  algorithms. 
At  each  level  of  refinement  and  for  each  problem,  all  calculations  were  completed  on  the  same 
computational  meshes.  Furthermore,  identical  convergence  criteria  were  applied  to  each  of  the 
algorithms. 

These  criteria  include  root  mean  square  (RMS)  measures  of  the  global  residual  of  all  conserved 
variables  as  well  as  integrated  surface  pressure. 


•  Global  residual: 


^--urhrMm 


where  k  denotes  the  kth  equation  on  an  IL  X  JL  computational  mesh  and  Rk  is  the  residual: 

Rk~  dt  ~  dx  dy 


For  plotting  purposes,  these  values  are  normalized  by  the  value  of  ||GM£.||  obtained  after 
the  first  iteration. 


•  RMS  Surface  pressure: 

IL 

Both  residuals  are  plotted  against  the  characteristic  time  Tc  =  L/V «>  where  L  is  a  characteristic 
length  scale  (chosen  as  the  cylinder  diameter)  and  is  the  freestream  velocity  of  the  incoming 
flow  upstream  of  the  impinging  shock.  Tc  provides  a  uniform  global  time  scale  for  comparison  of 
all  methods. 
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4.  Results 


4.0.1  Case  A  -  Blunt  Body 

Figure  2  displays  the  convergence  behavior  of  the  seven  schemes  for  Mach  8.03  flow  past  a  2-D 
cylinder.  The  calculations  progressed  until  the  residuals  dropped  five  orders  of  magnitude.  In 
this  calculation,  as  throughout  this  research,  no  emphasis  was  placed  on  efficiency.  Clearly,  by 
adjusting  some  parameter,  convergence  of  any  of  these  schemes  could  be  accelerated.  This  figure 
is  intended  only  to  demonstrate  that  all  schemes  converge  for  this  problem  and  provides  a  basis 
for  the  investigations  which  follow. 

Figure  3  exhibits  the  RMS  surface  pressure  for  each  scheme.  With  this  criterion,  all  schemes 
converge  within  one  characteristic  time  and  the  discrepancy  in  asymptotic  levels  may  be  ascribed 
to  slight  differences  in  normalization  and  are  considered  unimportant. 

Figure  4  shows  both  the  40  X  40  mesh  and  Mach  contours  for  each  of  the  seven  schemes. 
Mach  contours  are  drawn  in  increments  of  0.25  and  the  1.0  contour  (sonic  line)  is  annotated  in 
each  plot.  In  all  cases,  the  sonic  line  displays  the  highly  curved  profile  associated  with  2-0  flows 
above  Mach  2.  Additionally,  the  stand-off  distance  in  all  cases  agrees  with  the  value  published  in 
classical  texts  [42]  to  within  a  grid  point  (measured  from  the  “first  shock  point”  away  from  the 
body). 

While  the  contour  plots  give  some  indication  of  the  scheme’s  ability  to  capture  strong  shocks, 
a  more  precise  examination  stems  from  line  plots  of  Mach  number  along  specific  radial  mesh  lines. 
Figure  5  contains  such  a  plot  along  the  stagnation  streamline  while  Figure  6  contains  the  same 
plot  for  a  ray  inclined  36°  degrees  from  the  vertical.  These  plots  clearly  demonstrate  the  superior 
shock  resolution  of  the  upwind  schemes  as  compared  to  the  centered  methods.  Moreover,  they 
also  show  the  nonosdllatory  nature  of  the  flux  limited  upwind  schemes. 

Figures  7  and  8  show  the  same  line  plot  for  a  Mach  16.33  calculation  on  the  same  grid.  While 
the  central  difference  schemes  each  require  at  least  one  additional  grid  point  in  the  shock,  the 
upwind  schemes  still  resolve  the  shock  with  only  two  intermediate  cells.  Of  particular  note  are 


the  crisp  shocks  afforded  by  the  schemes  of  Roe,  van  Leer  and  Upwind  TVD.  The  relatively  poor 
performance  of  the  Ni  scheme  off- centerline  at  this  Mach  number  (Fig.  8)  may  be  ascribed  to  the 
use  of  the  diagonals  in  the  difference  stencil  -  in  contrast  to  the  other  schemes  which  excludes 
these  nodes  (e.g.,  [40]). 

Figure  9  displays  the  stagnation  enthalpy  losses  along  the  stagnation  streamline.  Following 
the  upwind  schemes  (3  through  7)  along  this  streamline,  perturbations  appear  within  the  discrete 
shock  and  all  curves  exit  with  slightly  different  (negative)  values.  Through  the  shock,  all  the 
flux  limited  schemes  reduce  their  approximation  of  the  governing  equations  to  first  order  as  the 
limiter  “clips  off”  the  solutions  at  extrema.  This  is  responsible  for  the  stagnation  enthalpy  error 
immediately  following  the  shock.  It  is  worth  noting  that  this  post-shock  error  decreases  for  all 
schemes  with  grid  refinement.  After  the  shock,  the  methods  again  become  second  order,  but 
despite  this,  only  the  van  Leer  and  Roe  schemes  remain  reasonably  constant  on  either  grid.  Of 
the  other  upwind  schemes,  only  the  Steger- Wanning  algorithm  is  of  particular  note  as  its  extreme 
diffusivity  creates  an  error  of  over  10%  on  both  grids.  In  evaluating  any  particular  method,  both 
its  overall  and  post-shock  behavior  must  be  considered.  Of  the  current  implementations,  the  van 
Leer  algorithm  shows  the  least  discrepancy  along  the  stagnation  streamline.  The  Jameson  central 
difference  scheme  (scheme  1)  exhibits  a  large  stagnation  enthalpy  drop  across  the  shock  due  to 
the  large  second  order  damping  required  to  capture  the  strong  Mach  8.03  normal  shock.  In  the 
subsonic  post-shock  region,  however,  considerable  improvement  is  observed  with  grid  refinement 
as  the  method  demonstrates  its  well  known  ability  to  solve  subsonic  flows  [43]. 

4.0.2  Case  B  -  Type  IV  Interference  Flow 

The  next  set  of  figures  presents  results  for  the  shock  interaction  test  case  at  Mach  8.03  on  a 
60  x  40  grid.  Figure  1  presents  details  of  the  geometry  and  parameters  specifying  this  problem. 
Figures  10  and  11  display  the  global  residual  and  surface  pressure  convergence  histories  for  all  the 
methods.  These  plots  display  converged  results  only  for  schemes  3,  5,  and  6  and  while  the  van 
Leer  and  Upwind  TVD  methods  exhibit  declining  residuals,  convergence  cannot  be  conclusively 
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claimed.  However,  the  surface  pressure  (Figure  11)  indicates  a  solution  converged  to  plotting 
accuracy.  Both  central  difference  schemes  appear  to  level  off  after  dropping  about  two  orders 
of  magnitude  and  show  small,  high  frequency  perturbations  around  a  constant  value.  However, 
these  oscillations  are  small  enough  that  no  changes  in  the  flow  field  are  observed.  These  schemes 
required  a  second  difference  smoothing  coefficient  about  1.5  times  that  needed  in  the  blunt  body 
calculation,  and  could  only  be  run  at  a  CFL  number  more  than  10  times  smaller  than  in  the 
noninterfering  case.  Notice  that  since  some  experimental  and  computational  evidence  exists  which 
suggests  an  unsteady  solution  all  schemes  were  run  with  a  constant  global  timestep  [23,  30,  25]. 

Figure  12  compares  contours  of  constant  Mach  number  for  the  schemes  under  investigation. 
These  contours  have  an  increment  of  0.25  and  are  labeled  in  the  last  frame.  This  figure  also 
contains  a  plot  of  the  coarse  computational  mesh  (which  was  not  optimized  excessively  for  this 
problem).  Proper  resolution  of  this  interaction  demands  both  accurate  representation  of  strong 
shocks  and  shear  layers  throughout  the  interference  pattern  within  the  distorted  bow  shock. 

As  anticipated,  the  additional  dissipation  required  to  curtail  oscillations  in  the  centered 
scheme’s  convergence  history,  considerably  smears  both  the  impinging  shock  and  the  distorted 
bow  shock.  All  the  upwind  methods  resolve  these  two  features  with  roughly  the  same  resolution. 
The  dissipative  nature  of  the  upwind  schemes  is  most  evident  in  the  Steger- Warming  splitting. 
Here,  the  internal  structure  of  the  field  is  smeared  even  more  than  with  the  central  schemes, 
and  the  contours  are  uniformly  smooth  throughout  the  field.  Of  the  four  schemes  (4  through  7) 
which  capture  the  embedded  shock  and  terminating  normal  shock,  the  van  Leer,  Roe  and  Upwind 
TVD  methods  perform  perceptibly  better  than  the  Symmetric  TVD  scheme.  All  four  algorithms 
predict  supersonic  jets  of  comparable  width. 

The  final  figure  for  the  coarse  grid  Type  IV  calculation  presents  normalized  surface  pressure  as 
a  function  of  9  (see  Fig.  1).  In  Figure  14,  the  vertical  axis  is  normalized  by  the  stagnation  pressure 
behind  a  Mach  8.03  normal  shock  (to  show  pressure  amplification).  Here,  results  from  the  seven 
calculations  are  compared  with  results  from  the  experiment  by  Wieting  and  Holden  [29].  Away 
from  the  peak,  all  schemes  agree  reasonably  well  with  the  available  data.  Each  of  the  schemes 
places  the  peak  at  roughly  the  same  angular  location,  but  are  also  shifted  about  3°  from  the 
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experimentally  observed  peak.  Since  the  angular  location  of  the  peak  is  a  strong  function  of  the 
impinging  shock  location,  this  small  discrepancy  is  not  considered  important  for  the  purposes 
of  the  present  investigation.  The  Steger- Warming  method  underpredicts  the  maximum  pressure 
augmentation  considerably.  Supporting  the  similarities  shown  in  the  contour  plots  of  the  previous 
figure,  the  Roe  flux  difference  splitting  and  the  Upwind  TVD  schemes  remain  nearly  inseparable 
on  this  plot,  achieving  a  value  somewhat  higher  than  the  Symmetric  TVD  scheme.  The  van  Leer 
splitting  reaches  a  slightly  higher  peak  but  is  shifted  roughly  2°  to  the  right  of  the  others.  Despite 
the  apparent  lack  of  structure  in  the  contour  plots  for  the  central  difference  schemes,  both  predict 
very  reasonable  values  for  the  peak  over  pressure. 

In  an  effort  to  bring  schemes  1,  2,  4  and  7  closer  to  convergence,  the  solutions  were  marched 
out  several  more  characteristic  times.  With  this,  both  the  van  Leer  and  Upwind  TVD  algorithms  > 

converged  fully.  However,  no  reasonable  combination  of  second  and  fourth  difference  dissipation 
coefficients  could  be  found  to  fully  lead  the  central  difference  schemes  to  convergence.  Over  the 
course  of  this  search,  the  peak  surface  pressure  varied  from  6  to  8.5  and  the  “loop”  structure  (ev¬ 
ident  in  the  solution  of  schemes  4  through  7  of  Figure  12)  occasionally  became  evident.  However, 
at  damping  coefficients  lower  than  those  reported  during  the  discussion  of  Figure  12,  the  high 
frequency  oscillations  shown  in  the  convergence  plots  increased  in  amplitude.  Figure  13  shows 
the  surface  pressure  history  of  the  Lax-Wendroff  scheme  with  a  second  difference  smoothing  coef¬ 
ficient  equal  to  that  of  the  blunt  body  case.  This  figure  shows  a  small  high  frequency  oscillation 
superimposed  on  a  large  periodicity  which  spans  about  10  characteristic  times.  This  plot  also 
shows  the  structure  of  the  flow  field  at  two  different  locations  in  the  pressure  cycle.  Note  that 
this  period  corresponds  to  a  physical  frequency  of  2200  Hz  which  is  slightly  lower  than  the  range 
(3  kHz- 10  kHz)  reported  by  Holden  et  al.  [23]. 

The  fact  that  this  unsteady  behavior  increased  as  diffusion  was  removed  from  the  solutions  led 
to  questions  about  the  diffusion  present  in  the  discrete  solutions  of  the  upwind  schemes.  Returning 
briefly  to  Figure  12  permits  a  crude  estimation  of  the  relative  diffusivity  in  the  upwinded  discrete 
solutions.  The  smooth,  rounded  contours  of  the  Steger- Warming  solution  are  indicative  of  excess 
dissipation  which  makes  itself  evident  in  the  underprediction  of  the  pressure  peak  in  Figure  14. 
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The  contours  in  the  Upwind  TVD  and  Roe  flux  difference  splitting  schemes  exhibit  nearly  identical 
smoothness  which  correlates  well  with  the  similarity  in  their  pressure  peaks.  The  van  Leer  result, 
however,  demonstrates  that  this  value  cannot  be  exclusively  linked  to  the  diffusion  in  the  field. 
Despite  the  fact  that  its  subsonic  contours  are  somewhat  smoother  than  schemes  5  or  7,  the 
pressure  peak  is  the  highest  of  all  the  upwind  algorithms. 

In  an  effort  to  make  a  more  quantitative  estimate  of  the  diffusivity  of  the  various  upwind 
methods  in  this  complex  flow,  the  truncation  error  in  the  solution  was  reduced  by  refining  the 
mesh.  Figures  15  and  16  contain  the  global  residual  and  surface  pressure  convergence  histories 
for  one  of  the  central  methods  (Scheme  1)  and  all  the  upwind  algorithms  on  a  120  x  80  mesh. 
On  this  finer  mesh,  only  the  convergence  behavior  of  the  two  flux  vector  split  schemes  remains 
unaltered.  Nevertheless,  the  surface  pressure  history  of  the  Symmetric  TVD  scheme  appears  to 
be  quite  steady.  Again,  all  these  calculations  maintained  a  constant  global  timestep  in  an  attempt 
at  time  accuracy.  Moreover,  marching  the  solution  out  several  characteristic  times  did  not  alter 
these  trends  significantly. 

Since  the  two  implicitly  integrated  flux  vector  split  schemes  (3  and  4)  were  the  only  ones  to 
converge  on  this  grid,  the  effect  of  the  time  integration  scheme  was  investigated  in  further  detail. 
The  van  Leer  splitting  was  reapplied  to  the  problem  using  a  second  order  time-accurate  explicit 
predictor-corrector  algorithm  on  both  meshes.  Both  convergence  criteria  were  unaffected  on  the 
coarse  mesh  (Fig.  17).  On  the  medium  grid,  however,  the  pressure  residual  remained  unchanged 
but  the  global  residual  with  the  explicit  integration  dropped  only  two  orders  of  magnitude  and 
then  remained  generally  steady  -  albeit  with  small  amplitude  high  frequency  oscillations.  No 
difference  in  surface  pressure  or  flow  structure  could  observed  at  these  levels. 

Figure  18  displays  the  contour  plots  of  constant  Mach  number  on  the  120x80  mesh.  Of  course, 
the  detailed  structure  in  the  plots  for  the  oscillatory  behaved  solution  is  a  function  of  the  current 
characteristic  time.  Nevertheless,  all  the  plots  presented  are  representative  pictures  of  the  solution 
provided  by  each  scheme.  All  the  methods  show  increased  resolution  as  compared  to  the  coarse 
grid  plots  of  Figure  12.  On  this  medium  grid,  both  the  Jameson  and  Steger- Warming  schemes 
show  the  loop  structure  demonstrated  by  the  other  upwind  schemes  on  the  coarse  grid.  Despite 
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the  accurate  shocks  evident  in  the  Steger- Warming  calculation,  the  internal  structure  is  less 
apparent  than  in  the  central  difference  solution.  Again,  methods  4  through  7  demonstrate  superior 
resolution  in  these  plots.  As  on  the  coarse  mesh,  the  van  Leer  and  Symmetric  TVD  methods  give 
very  similar  looking  results.  After  passing  through  the  terminating  normal  shock,  most  of  the 
flow  accelerates  and  expands  upward  around  the  body.  This  flow  forms  a  supersonic/subsonic 
shear  layer  with  the  nearly  stagnated  flow  behind  the  bow  shock.  Both  the  Roe  and  Upwind  TVD 
methods  appear  to  show  a  developing  shock  reflection  system  within  this  shear  layer. 

Figure  19  presents  a  detailed  view  of  a  typical  result  (Upwind  TVD)  through  velocity  vectors. 
Figure  20  shows  a  few  representative  particle  traces  (PLOT3D)  superimposed  on  selected  Mach 
contours.  The  streamlines  detail  the  flow’s  progress  through  the  interference  and  stagnation 
regions  of  the  flow  field. 

Figure  21  contains  a  plot  of  normalized  surface  pressure  for  the  methods  on  the  medium 
grid.  Note  again  that  the  unsteady  result  for  the  Roe  flux  difference  splitting  and  Upwind  TVD 
methods  are  shown  as  “snapshots”  taken  at  the  same  characteristic  time  as  the  contours  in 
Figure  18.  Qualitatively,  all  the  peaks  on  this  plot  are  higher  than  on  the  coarse  grid. 

As  a  further  investigation  of  the  oscillatory  behavior  demonstrated  by  the  Upwind  TVD  and 
Roe  flux  difference  split  schemes,  Figures  22  and  23  present  the  solution  of  the  Upwind  TVD 
algorithm  at  extrema  in  the  solution  cycle.  Both  low  and  high  frequency  oscillations  appear 
in  the  surface  pressure  convergence  behavior.  The  convergence  history  is  reproduced  on  both 
figures  for  the  solution  between  90  and  120  characteristic  times.  While  there  is  little  evidence  of 
the  low  frequency  oscillations  in  the  Mach  and  surface  pressure  distributions,  the  high  frequency 
oscillation  has  dramatic  effects  on  both.  The  unsteadiness  manifests  itself  through  changes  in  the 
angle  of  the  terminating  normal  shock,  which  in  turn  influences  the  angle  at  which  the  jet  impinges 
on  the  surface.  This  sets  up  different  oblique  shock  reflection  patterns  within  the  supersonic  layer 
bounded  by  the  sub  /supersonic  shear  layer  and  the  body.  Looking  more  closely  at  Figures  22 
and  23,  the  higher  peaks  in  the  surface  pressure  plot  correspond  to  the  pattern  where  the  jet 
terminating  normal  shock  is  more  closely  aligned  with  the  body  surface.  Returning  to  the  solution 
with  the  van  Leer  scheme  on  this  grid  (Figures  18  and  21),  both  the  terminating  shock  angle  and 
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the  surface  pressure  distribution  with  this  method  correspond  very  closely  to  the  Upwind  TVD 
solution  at  the  maximum  point  in  the  oscillation  of  the  RMS  surface  pressure  distribution. 

Although  the  van  Leer  and  Symmetric  TVD  schemes  provided  reasonably  accurate  surface 
pressure  distributions  (Figure  21),  they  fail  to  resolve  the  reflected  shock  structure  within  the 
supersonic  region  bounded  by  the  upper  body  and  the  sonic  line  (Figure  18)  seen  in  the  discrete 
solutions  of  Roe  and  Upwind  TVD,  Moreover,  neither  method  predicted  an  unsteady  flow  solution. 
With  the  intent  of  resolving  more  structure  within  the  supersonic  layer  and  to  further  investigate 
the  convergence  behavior  of  the  van  Leer  and  Symmetric  TVD  schemes,  the  truncation  error  was 
reduced  further  by  again  refining  the  mesh.  These  two  schemes  were  applied  to  a  mesh  constructed 
by  doubling  the  number  of  points  in  each  coordinate  direction  (resulting  in  a  240  x  160  mesh). 

Figures  24  and  25  present  convergence  behavior  and  Mach  contour  distributions  for  the  van 
Leer  and  Symmetric  TVD  discrete  solutions  on  the  fine  grid.  In  Figure  24  the  upper  frame 
chronicles  the  global  residual  variation  for  50  characteristic  i  imes.  As  anticipated,  both  schemes 
now  predict  unsteady  behavior  with  high  frequency  oscillations  as  seen  earlier  in  Roe  and  Upwind 
TVD  solutions  on  the  medium  grid  (Figures  15  and  16).  The  lower  frame  records  RMS  surface 
pressure  over  the  same  period  of  time.  The  smaller  amplitude  of  the  oscillations  of  the  van  Leer 
scheme  correspond  to  the  lower  magnitude  of  its  residual  in  the  upper  frame  of  this  figure.  In 
Figure  25,  both  these  schemes  now  capture  the  reflected  shock  structure  within  the  supersonic 
layer  flowing  up  over  the  body  surface  as  seen  on  the  medium  grid  with  the  Roe  and  Upwind 
TVD  methods.  With  reduced  truncation  error  present  on  the  fine  grid,  both  solutions  are  now 
unsteady  and,  as  a  result,  slight  differences  may  be  expected  between  the  two  plots  despite  higher 
order  time  accuracy. 

Figure  26  presents  a  schematic  of  the  flow  structures  evidenced  by  the  discrete  solutions  in 
Figures  18,  22  and  25.  This  schematic  identifies  the  supersonic  flow  region  bounced  by  the  body 
and  the  sub/supersonic  shear  layer  which  forms  around  the  upper  surface  of  the  body.  The 
picture  suggests  a  possible  feedback  mechanism  which  may  be  responsible  for  the  unsteadiness  in 
the  oblique  shock  reflection  patterns  shown  in  Figure  22.  With  the  terminating  normal  shock  at 
some  initial  angle,  an  associated  oblique  shock  pattern  will  establish  itself  between  the  shear  layer 
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and  the  body.  As  these  oblique  shocks  form,  they  disturb  this  shear  layer,  communicating  these 
disturbances  to  the  subsonic  region.  Within  the  subsonic  stagnation  region,  these  disturbances 
are  free  to  propagate  in  all  directions  including  back  toward  the  upper  shear  layer  (see  Fig.  1).  As 
these  waves  distort  the  interaction  region,  they  may  change  the  angle  of  the  terminating  normal 
shock  setting  up  a  new  oblique  shock  reflection  system  which  again  distorts  the  subsonic  shear 
layer  thus  reinitiating  the  cycle.  While  this  hypothesis  seems  reasonable  and  physical,  it  does  not 
satisfactorily  explain  the  difference  in  frequency  between  the  computed  and  experimental  results. 
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5.  Conclusion 


A  comparative  study  has  been  performed  for  a  variety  of  upwind  schemes  applied  to  a  problem 
of  realistic  complexity.  The  schemes  which  achieve  second  order  accuracy  through  the  MUSCL 
approach  included  the  Steger  Warming  and  van  Leer  flux  vector  split  schemes  and  the  Roe  flux 
difference  split  scheme.  Additionally,  this  comparison  investigated  two  non-MUSCL  algorithms 
based  on  the  “Symmetric”  and  “Upwind”  TVD  methods  of  Yee  and  Harten  and  Yee. 

The  shock  capturing  ability  of  the  upwind  schemes  was  demonstrated  for  grid  aligned  and 
nonaligned  shocks  at  two  different  Mach  numbers.  This  investigation  clearly  showed  the  superior 
shock  resolution  afforded  by  these  procedures.  Moreover,  this  ability  was  found  to  be  independent 
of  Mach  number,  which  is  not  the  case  for  central  schemes. 

The  Type  IV  shock  interference  flow  pointed  out  differences  among  the  upwind  procedures.  In 
particular,  the  extreme  diffusivity  of  the  Steger- Warming  method  -  which  was  not  immediately 
evident  in  the  blunt  body  solutions  -  became  apparent  and  unacceptably  corrupted  the  solution. 
The  other  four  upwind  procedures  all  captured  considerable  flow  field  detail  with  a  slight  advantage 
falling  to  the  Upwind  TVD  and  Roe’s  algorithm. 

A  comparison  of  flow  structure  on  a  successively  refined  sequence  of  three  meshes  indicated 
the  least  diffusion  in  the  Upwind  TVD  and  Roe  flux  difference  split  discrete  solutions.  These 
two  schemes  already  showed  oscillatory  behavior  on  a  120  x  80  grid.  The  Symmetric  TVD  and 
van  Leer  algorithms  also  led  to  oscillatory  flows  for  this  problem  but  required  an  even  finer  mesh 
(240  x  160)  to  demonstrate  this  behavior. 

Despite  the  apparent  lack  of  structure  within  the  flow  solution  caused  by  the  inability  of  the 
central  difference  algorithms  to  resolve  features  on  the  order  of  the  mesh  scale,  both  these  schemes 
provided  reasonable  surface  pressure  distributions.  Moreover,  these  procedures  demonstrated 
neutral  convergence  behavior  on  the  coarse  grid.  Although  such  behavior  became  apparent  in 
some  upwind  procedures  on  finer  meshes,  it  is  not  clear  that  the  same  mechanism  is  responsible 
for  the  unsteadiness  observed  in  the  cent^ ;  and  upwind  solutions. 
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Nomenclature 

a  local  sound  speed 

A  flux  Jacobian  matrix  (f$) 

f\  van  Leer  mass  flux  function  for  subsonic  flow 

g  limiter  function  for  “Upwind”  TVD  scheme 

P,  G  flux  density  vectors 

G.R.  global  residual 

IL ,  JL  number  of  mesh  points  in  each  coordinate  direction 
L  .  characteristic  length 

M  Mach  number 

p  pressure 

Q  numerical  flux  function  in  TVD  schemes 

R  gas  constant,  right  eigenvector  matrix,  residual 

t  nondimensional  time 

T  temperature 

Tc  characteristic  time 

u,  v  Cartesian  velocity  components 

U  vector  of  conserved  variables 

V  magnitude  of  total  velocity  vector 

x,  y  two-dimensional  Cartesian  coordinates 

a  difference  in  characteristic  variables 

6  small  constant  6  <  1 

A  difference  between  left  and  right  states 

A  limiter  function  for  MUSCL  type  schemes 
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7 

ratio  of  specific  heats 

A 

eigenvalue 

A 

eigenvalue  matrix 

<t> 

element  of  $ 

$ 

limited  flux  function  in  TVD  schemes 

9 

circumferential  angle 

P 

density 

V> 

function  enforcing  entropy  condition 

Subscript 

(  M 

face  to  the  right  of  tth  node 

(  )x 

x  direction 

(  )oo 

freestream  value 

Superscript 

(  )' 

/th  component  of  vector 

(  t 

state  to  left  of  current  face 

(  f 

state  to  right  of  current  face 

(  )+,_ 

positive  and  negative  components 

C  ) 

Roe  averaged  state 
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Distorted  Bow  Shock 


Region  1  Region  2 


Figure  1:  Schematic  of  Type  IV  shock-shock  interaction  of  Edney  [28].  Angies  of  shocks  and 
shear  layers  calculated  from  oblique  shock  relations  and  hodograph  methods  of  [44]. 
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Figure  4:  Mach  contours  for  blunt  body  on  a  40  x  40  grid  at  Mach  8.03. 
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Figure  5:  Mach  number  vs.  grid  point  for  a  blunt  body  on  a  40  x  40  grid  at  Mach  8.03  along 
the  centerline. 


Coarse  (20  x  20) 


009 


pressure  convergence  behavior  for  Type  IV  interaction  on  coarse  (60  X 
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Figure  12:  Mach  contour  plots  for  Type  IV  interaction  on  coarse  (60  x  40)  grid,  Moo  =  8.03 
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Figure  13:  Mach  contours  at  various 

with  dissipation  reduced  to  blunt 
grid,  Moo  =  8.03). 


Figure  14:  Surface  pressure  comparison  for  Type  IV  interaction  on  coarse  (60  X  40)  grid  M< 
8.03. 
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RMS  surface  pressure  convergence  behavior  for  Type  IV  interaction  on  medium 
(120  x  80)  mesh,  -  8.03. 


Symmetric  TVD  Upwind  TVD 

Figure  18:  Mach  contour  plots  for  Type  IV  interaction  on  medium  (120  x  80)  grid,  =  8. 
(continued) 


Fip.re  19:  Velocity  vector,  of  typical  solution  of  Type  IV  interaction  on  a  medium  (120  x  80) 
grid  (Upwind  TVD),  Af*,  =  8.03. 
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Figure  21:  Surface  pressure  comparison  for  Type  IV  interaction  on  medium  (120  x  80)  grid, 


Figure  22:  Mach  contours  at  various  extrema  in  convergence  cycle  for  Upwind  TVD  scheme 
medium  (120  x  80)  grid,  Moo  —  8.03. 


various  extrema  in  convergence  cycle  for  Upwind  TVD  scheme 


RMS  Surface  Pressure  Residual 


Characteristic  time 


Figure  24:  Convergence  behavior  for  Type  IV  interaction  on  fine  (240  x  160)  grid,  Moo  =  8J 
Upper:  Residual,  Lower:  RMS  surface  pressure. 
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